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ABSTRACT 

An efficient implementation of the generalized method of cells micromechanics model is presented 
that allows analysis of periodic unidirectional composites characterized by repeating unit cells containing 
thousands of subcells. The original formulation, given in terms of Hill’s strain concentration matrices that 
relate average subcell strains to the macroscopic strains, is reformulated in terms of the interfacial subcell 
tractions as the basic unknowns. This is accomplished by expressing the displacement continuity equa- 
tions in terms of the stresses and then imposing the traction continuity conditions directly. The result is a 
mixed formulation wherein the unknown interfacial subcell traction components are related to the 
macroscopic strain components. Because the stress field throughout the repeating unit cell is piece-wise 
uniform, the imposition of traction continuity conditions directly in the displacement continuity equa- 
tions, expressed in terms of stresses, substantially reduces the number of unknown subcell traction (and 
stress) components, and thus the size of the system of equations that must be solved. Further reduction in 
the size of the system of continuity equations is obtained by separating the normal and shear traction 
equations in those instances where the individual subcells are, at most, orthotropic. The reformulated ver- 
sion facilitates detailed analysis of the impact of the fiber cross-section geometry and arrangement on the 
response of multi-phased unidirectional composites with and without evolving damage. Comparison of 
execution times obtained with the original and reformulated versions of the generalized method of cells 
demonstrates the new version’s efficiency. 


INTRODUCTION 

The two-dimensional generalized method of cells is a micromechanics model developed by Paley 
and Aboudi (1992) and Aboudi (1993) for predicting the response of unidirectional metal matrix compo- 
sites with periodic microstructures. A continuously reinforced composite is modeled as a doubly periodic 
assemblage of fibers embedded in a matrix phase, Fig. 1. The periodic character of the assemblage allows 
one to identify a repeating unit cell that is the building block for the entire composite. The properties of 
this repeating unit cell are thus representative of the properties of the entire assemblage. The repeating 
unit cell consists of NpxN y subcells. Each of these subcells is assumed to be occupied, in general, by a 
material that exhibits inelastic behavior. The subcell material’s inelastic behavior can be modeled by a 
variety of constitutive theories, including linear viscoelasticity, classical incremental plasticity, or unified 
viscoplasticity theories. Thus the repeating unit cell consists of N pxNy different inelastic materials, i.e., it 
represents a multi-phased, inelastic composite. 

The generalized method of cells allows approximate micromechanical analysis of fairly complicated 
periodic arrays, including: 
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• thermomechanical response of multi-phased, metal matrix composites 

• modeling of variable fiber shapes 

• analysis of different fiber arrays 

• modeling of porosities and damage 

• modeling of interfacial regions around inclusions, including interfacial degradation 


The effective or average stress-strain equations for the composite are constructed through an 
approximate deformation analysis in each subcell of the repeating unit cell based on a linear representa- 
tion of the displacement field in terms of local subcell coordinates. Using this approach, the governing 
field equations within each subcell are satisfied identically in a volumetric sense, while the traction and 
displacement continuity conditions between adjacent subcells are imposed in an average sense together 
with an homogenization condition that ensures that the response of a given repeating unit cell is indistin- 
guishable from that of its neighbors. The above analysis establishes the so-called Hill’s strain concentra- 
tion matrix relationships between the average subcell strains and the imposed average or composite 
strains, that are used in the construction of the composite stress-strain equations (Hill, 1963). The size of 
the strain concentration matrix is 6N$Ny x 6NpNy since there are six unknown subcell strains in each of 
the NpN y subcells. This matrix must be inverted once if an isothermal mechanical analysis is conducted, 
and many times if a nonisothermal thermomechanical analysis is conducted with temperature-dependent 
constitutive properties in the individual subcells. Consequently, this imposes limits on the level of 
discretization of the repeating unit cell. Practically, the analysis of repeating unit cells larger than 10x10 
may become computationally inefficient, particularly if the micromechanical analysis is part of a larger 
structural analysis. 

To enhance the computational efficiency of the generalized method of cells, Orozco (1997) took 
advantage of the sparse features of the strain concentration matrix in inverting it. The sparse implementa- 
tion of the generalized method of cells made possible the elastic analysis of periodic fiber arrays charac- 
terized by repeating unit cells containing thousands of subcells. Examples were presented for the calcula- 
tion of the effective properties of unidirectional composites containing up to lOOx 100 subcells in the 
repeating unit cell. 

Herein, another approach is presented for the efficient implementation of the generalized method of 
cells. This approach involves reformulating the displacement continuity equations in terms of the interfa- 
cial subcell traction components as the basic unknowns in place of the subcell strain components. This 
reformulation, together with the piece-wise uniform character of the stress field throughout the repeating 
unit cell and direct imposition of the interfacial traction continuity conditions between the individual sub- 
cells, substantially reduces the size of the matrix that must be inverted to express the unknown interfacial 
subcell tractions (and thus subcell stress components) in terms of the macroscopic strains. This matrix 
consists of diagonal submatrices and fully-populated off-diagonal submatrices, thereby lending itself to 
further computational efficiency enhancement through Orozco’s sparse implementation approach. 
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STANDARD FORMULATION OF THE GENERALIZED METHOD OF CELLS 

Since the reformulation of the generalized method of cells requires the use of the original equations, 
we begin by providing a brief outline of the original formulation. We proceed to outline the original 
micromechanical analysis by first defining the relationships between the average composite and the aver- 
age subcell stress and strain quantities. 

The volume-averaged subcell stress is defined in the usual way, 

+/lp/2 +ly/2 

g» = J_ J j otw.ef 4* (1) 

v ¥>y 2 -1/2. 

(where Vp T = hpl y ). Thus the average composite stress a is obtained from a weighted sum of the subcell 
average stresses taken over all the subcells. 


a = 


1 v v/, / s (pY) 

ttZEVy 0 

m p=ir=i 


( 2 ) 


where h and l are the dimensions of the repeating unit cell, Fig. 1. Similarly, the volume-averaged subcell 

. _ (Pt) . . , 

strain e is given by, 


+/y /2 

_(Py) _ _L J j e (Py )dx2&f (3) 

^PY -ftp/2 - 1/2 

and the average composite strain e is obtained from a sweighted sum of the subcell average strains taken 
over all the subcells, 


£ = 


] N f N i (Py) 

m P=1t=1 


(4) 


The relationship between average subcell stresses and strains is obtained by volume averaging the 
constitutive equations used to describe the material behavior in the subcell ((3y), 


a m = C (Py) [e^ -t 9 * - a (pY) AT ] 


(5) 


where I p(Pt) is the average plastic subcell strain, C (pY) is the elastic stiffness matrix, a (Pr) are the coeffi- 
cients of thermal expansion, and AT is the temperature deviation. The above description of the material 

. _p(py) . , 

behavior is sufficiently general to admit any inelastic constitutive model for the plastic strain e in the 
subcell (py). Herein, the classical incremental plasticity theory, reformulated in terms of strains (Mendel- 
son, 1983; Williams and Pindera, 1997) is used to model the inelastic effects. 
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The effective or average stress-strain equations for the composite are constructed from the definition 
of the volume-averaged composite stresses, eqn (2), by first expressing the volume-averaged subcell 
stresses, eqn (1), in terms of the volume-averaged subcell strains using eqn (5). Then, the volume- 
averaged subcell strains are expressed in terms of the composite strains through the use of Hill’s strain 
concentration matrix relations obtained from the approximate deformation analysis of the individual sub- 
cells. This produces the effective stress-strain equations for the composite in the form: 


o = C* (e-f-aAT) (6) 

where f? is the average composite plastic strain, C* is the effective elastic stiffness matrix, and a* is the 
effective thermal expansion coefficient vector. 

The relationship between the volume-averaged subcell and composite strains, necessary to generate 
eqn (6), is obtained through an approximate analysis of the deformation fields in each subcell of the 
repeating unit cell. Towards this end, the displacement field in the individual subcells modeling a con- 
tinuously reinforced composite is approximated in terms of a linear expansion in the local coordinates 
xf, xf centered at the mid-point of a given subcell, 

uffi = co^> + xf + xf yjw (7) 


where i = 1, 2, or 3; (3 = 1, ..., jVp; y= 1, ..., N y \ and co^ are the displacement components at the center 
of each subcell ((3y), defined at discrete points in the x 2 -jc 3 plane but treated as continuous functions of 
the axial coordinate x i aligned with the fiber direction. Using the local strain-displacement equations for 
each subcell, one can easily show that the micro variables <j>^ and are related to the subcell strains 
e^. It follows then, that the subcell strains, and therefore stresses, are piece- wise uniform throughout the 
repeating unit cell. Thus we obtain the following relations between the microvariables <|>P Y) and and 

_(Py) 

the volume-averaged subcell strains £ /; , 



— (Py) 
£ 12 



dx\ 


OX I 

2 ox \ 2 


(8a) 

(8b) 


Since the subcell stresses and strains are piece-wise uniform throughout the repeating unit cell, the 
governing field equations are satisfied identically. The piece-wise uniform character of the stress and 
strain fields, however, requires that the interfacial traction and displacement continuity between the indi- 
vidual subcells within the repeating unit cell, as well as between the given repeating unit cell and the sur- 
rounding cells, be imposed in an average sense. Imposing these displacement continuity conditions and 
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applying an homogenization condition that ensures that the given cell deforms in an indistinguishable 
manner from its neighbors, we obtain: 




P=l,...,Np 


where it can be shown that (Brayshaw, 1994) 


8<ojh> _ M 

d.Xj d.Xj 


(9a) 

(9b) 


(10) 


for all the subcells. Use of the above relations in the definitions for the average composite strains given in 
terms of the volume-averaged subcell strains, eqn (4), together with eqns (8) and (9) yields the relations. 


1 3 ( 0 , d(Oj 
£,J ~ 2 ( dxj + dx, ) 


(ID 


Thus the uniform axial deformation constraint and the interfacial displacement continuity conditions pro- 
vide the following relations between the volume-averaged subcell strains and the composite strains. 


_(Py) 

e n 


-Ell? P !>•••» 


N p ; y= 1,-, 


N y 


(12a) 


(By) 

E*p£22 =fce 22 , y= 1,..., Ny 

p=i 


(By) 

S>pe 12 =hl n , y= 1,..., N y 

p=i 


(12b) 

(12c) 


X/ y E 33 Y) =/E33, P=1 (Vp ( 12d ) 

Y= 1 


^Y (R v ) 

X/ 7 e,3 =(Ei3. P= 1,..., Np 

7-1 


(12e) 




"t (Pv) 

£ W y E23=W£ 23 

P=lY=l 


(12f) 
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This system of equations can be written in the matrix form, 

^g £ j = -^e (13) 

The matrix Ac contains terms that describe the internal geometry of the repeating unit cell. The vector 
z s = ( £* U \ J ? N * Ny) ) with dimension 6NpN y contains the six volume-averaged strain components for 
each of the N$N y subcells. The matrix J contains the overall cell dimensions. 

The displacement continuity conditions, eqns (13), provide 2(N$ + N y ) + N$N y + 1 equations 
involving the 6N$N y unknown subcell strains. The imposition of the interfacial traction continuity pro- 
vides the remaining 5N$N y - 2(Np + N y ) - 1 equations for the determination of the unknown subcell 
strains in terms of the macroscopic strains, 

P = P+1 when P < Np and P = 1 when P = N$, y= 1 N y (14a) 

Y = Y+1 when y <N y and y=\ when y =N y , P=l,.„, Np (14b) 

for i=l,2,3. These traction continuity equations can be expressed in matrix form in terms of volume- 
averaged total, plastic and thermal subcell strains upon use of the volume-averaged subcell constitutive 
relations, eqn (5), 

A M (e,-e?-e[) = 0 (15) 

, -P ,-P(ll) _p(AVV 7 ) -T _T(1I) _T(Ayv r ). . , . . ,. . 

where e s = ( e , ..., e ), e s = ( e , ..., e ), are the plastic strain and thermal strain vec- 
tors, respectively, for all the subcells, = a.^ AT, and the matrix Am contains the elastic constants of 

the materials in the individual subcells. 

Combining eqns (13) and (15), the system of equations relating the subcell strains to the uniformly 
applied composite strains, obtained from the imposition of the continuity conditions, can be written in 
compact matrix notation as follows, 

Az 5 = kl + b{€ s +z] ) (16) 

where 



The above equations can be solved in order to express the subcell strains in terms of the average strains: 
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( 18 ) 


e s = Ae + D(e^ + eJ) 

where A = A _1 K and D - A D. The matrix A, in fact, is the resulting Hill’s concentration matrix which 
relates the microstrains to macrostrains. It follows from eqn (18) that the average subcell strains can be 
expressed in terms of the macroscopic strains and average plastic and thermal subcell strains as follows: 

+ d 9 ) 

so that the average subcell stresses become: 

o (P7) = c^tA^e + D^{ t + e[ ) - +e' m )] (20) 

Using these expressions in the definitions for the macroscopic stresses given in terms of the average sub- 
cells stresses, eqn (2), we obtain the following constitutive equation for the composite response, 

a = - X ZHh c(H [A (PY) e + D (Py) ( e? + z , ) - (?*** + e^)] (21) 

M p=i-p=i 

Comparison of eqns (21) and (6) yields the following expressions for the macroscopic stiffness matrix C 
given in terms of the subcell stiffness matrices C (Py) and the subcell strain concentrations matrices A (Py) , 

= ZVy c( ^ a(N (22) 

hl p=t y=\ 

and for the macroscopic plastic and thermal strains e P and given in terms of the average subcell plastic 
and thermal strains e P<Pr> and g 7 ^ and the plastic strain concentration matrices D (Py) , 

t = -^(Cy^i^hpl y C^\D^ s (23) 
hl 

^ = -T7 (CV (24) 

ht P=lY=l 

Since the evolution of subcell plastic strains depends on the deformation history of the composite, these 
strains must be determined in an incremental fashion dependent on the particular form of the employed 
inelastic model for the matrix phase. A detailed presentation of the two-dimensional generalized method 
of cells analysis can be found in Paley and Aboudi (1992) and Aboudi (1993). 
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EFFICIENT REFORMULATION OF THE GENERALIZED METHOD OF CELLS 

The core of the computational effort in the original formulation of the generalized method of cells 
lies in evaluating the strain concentration matrices A and D in eqn (18), i.e., solving the 6N$N y x 6 NpN y 
system of equations given in eqn (16). The computational effort can be tremendously reduced by express- 
ing the displacement continuity relations, eqns ( 12 b)-( 12 f), in terms of subcell stresses using the strain- 
stress equations, 

= S (Py) + <x (P y) AT (25) 


and then imposing the interfacial traction continuity conditions directly in the reformulated interfacial 
displacement continuity equations. Further reduction can be achieved by separating the relationships 
between the subcell normal tractions and macroscopic normal strains from the relationships between the 
subcell shear tractions and macroscopic shear strains. When the subcells are, at most, orthotropic (as is 
the case here), this separation is possible because of the absence of shear-normal coupling owing to the 
use of a first order displacement expansion and the imposition of the interfacial continuity conditions in 
an average sense. The relations between subcell normal tractions and macroscopic normal strains are 
given first, followed by the corresponding shear relationships. 

Relations between subcell normal tractions and macroscopic normal strains 

We express eqns (12b) and (12d) in terms of the three subcell normal stresses using eqn (25). The 

(p-y) _ 

axial normal subcell stress o u is then expressed in terms of the common macroscopic strain e n and the 

_(p Y ) _(py) 

transverse normal subcell stresses O 22 and <*33 using the uniform axial deformation constraint, eqn 
(12a), and the first of the strain-stress equations (see the Appendix). Defining the common interfacial nor- 
mal traction 7$ for a fixed column of subcells (ly) (N$y) and for a fixed row of subcells 
((31) (pNy) dictated by the interfacial traction continuity conditions, 


— (ly) _(2y) 
^22 -<?22 


= y= 1,..., N y 


_(P1) _ _(p2) 
<7 33 _ o 33 


0^ = 7$, p=l,...,N|j 


(26) 

(27) 


the reformulated interfacial displacement continuity conditions become. 


A y T$ + =ht 2 2 -Cy€n +d y AT +p\ Y) , y= 1,..., N y 

p=l 


^,l y Bp y T$ +DpT$ =/e 33 -epEn +/pA 7 ' +p ^ , P= 1,..., Np 


(28) 

(29) 
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Equations (28) and (29) can be cast in the matrix form given below. 


A B 
B' D 


rk{ ok + j °k “{ k" + { /1 A7 ' + I A 


( 30 ) 


where A,B,B' and DaieN y x N y , N y x Np, Np x N y , and N p x Np matrices, respectively, whose struc- 
ture and elements are given below; T 2 = [7$ kk ] and 7 3 = [7$ Tjk] contain the N y and N p 
unknown normal tractions in the individual columns and rows of the unit cell, respectively; H and L are 
N y x 1 and Np x 1 vectors whose elements are the cell dimensions h and /, respectively; and the vectors 
c = [c c Ny ], d = [d\ d Ny ], e- [e i e Nf \, f- \f\ p \ - [M 0 »— » 

P2 = \P? , p-k ] contain elements also provided below. 


A, 

0 ■ 

• 0 


Mu 

^2^21 


0 

A 2 • 

• 0 

B = 

M,2 

h 2 B 2 j 


0 

0 • 

* 

1 


B\B ijv t 





/ 2 Si2 

i 

sT 

0Q 

►- 

J 


D, 

0 

• 0 

l,B u 

l 2 B 2 2 



0 

d 2 

• 0 

• 

• 


D = 

• 

- 

• 

^l^Npl 

^2^jVp2 



0 

0 • 



Vp s$?) 2 


), = 


Spsp 

sW> ’ 


Dp = £/ T (S^ 

t=i 


w ) 


C 


y 



d y 



-p(Py), 
-£22 ) 


N, 

ep = 

r=i 



V T rtfcf) 
fp = SZ^-SrOV 

Y=l 




(&YK (B> ✓‘^^(Py) sP(PDx 

- a 3 r ) » = -£33 ) 


(31a) 


(31b) 


(31c) 


Equation (30) forms the core of the computational effort in the reformulated generalized method of 
cells. The system of equations to be solved contains Np + N y unknowns in place of the 6NpN y unknowns 
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in the original version. Additional relations between subcell shear tractions and macroscopic shear strains 
do remain to be established to complete the reformulation. However, the computational effort involved in 
calculating the subcell shear tractions in terms of the macroscopic shear strains is minimal due to the 
absence of shear-normal and shear-shear coupling as illustrated next. Furthermore, we note that eqn (30) 
can be further reduced by expressing it either in terms of T 2 or 7' 3 . However, as will be illustrated in the 
next section, this reduction is offset by a greater number of matrix operations necessary to calculate the 
effective composite response, thus leading to increased efficiency only in certain circumstances. 

Relations between subcell shear tractions and macroscopic shear strains 

To obtain the axial shear relations, we express eqns (12c) and (12e) in terms of the subcell axial 
shear stresses using eqn (25). Owing to the absence of normal-shear and shear-shear coupling in orthotro- 
pic materials, these relations involve only the subcell shear stress and macroscopic shear strain quantities 
in the respective planes (see the Appendix). Defining the common interfacial shear traction T$ for a 
fixed column of subcells (ly) ,..., (N py) and for a fixed row of subcells (Pi) ,..., (P N y ) dictated by 

the interfacial traction continuity conditions, and utilizing the symmetry of the stress tensor, 

02 ? = 02? = ... = o ( 2 ? y) = 7$ =7$ , y= 1,..., N y (32) 

of? = of? = ... = = TfS= 7# , P = 1,-, Nfi (33) 

the reformulated interfacial displacement continuity conditions become, 

|(I/ip^)^=he 12 -Sh p ^? Y) , y= 1,..., N y (34) 

1 P=1 

^(1^)7$ = /e !3 - , P= I,-, tfp (35) 

1 1=1 7=1 

Thus the solution for the subcell interfacial shear tractions 7$ and is readily obtained in terms of the 
respective macroscopic shear strains and subcell plastic shear strains with minimum computational effort. 

To obtain the remaining transverse shear relation, we express eqn (12f) in terms of the subcell 

transverse shear stresses using eqn (25). Owing to the absence of normal-shear and shear-shear coupling 

in orthotropic materials, this relation involves only the subcell shear stress and macroscopic shear strain 
quantities in the transverse plane (see the Appendix). Defining the common interfacial shear traction T$ 
for a fixed column of subcells (ly) ,..., (Npy) and for a fixed row of subcells (Pi) ,..., (P N y ) dictated 
by the interfacial traction continuity conditions, 
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-dy) -(2y) 

<*23 ~ °23 


( 36 ) 


_(PD -(P2) 
<*32 “ <*32 



II 

N y 

= 5f’’ = 7f. 

II 

... 


(37) 


and utilizing the stress tensor’s symmetry, i.e., O 23 ~ ®i 2 » so that T’S? ~ 7^2 — T 23 f°r all combinations 
of P = l,...,iVp and y= 1,...,7V Y , the reformulated interfacial displacement continuity condition becomes, 

. Nf Ny Nf Ny _n(Bv) 

\(I,W y Sp)T 23 = hll 23 - (38) 

1 (3=ly=l P=ly=l 


Macroscopic constitutive equations 

The macroscopic constitutive equations for the composite are readily determined in the form of eqn 
(6) by first solving eqns (30), (34), (35), and (38) for the unknown interfacial subcell tractions (and thus 
subcell stresses), and then using the resulting expressions for the subcell stresses given in terms of the 
macroscopic strains and subcell plastic strains in the definitions for the macroscopic stresses, eqn (2). To 
determine the explicit expressions for the effective stiffness moduli, plastic and thermal strains in terms 
of the individual subcell moduli, plastic strains, thermal expansion coefficients and geometry based on 
the reformulated approach, we start with the solution of eqn (30). The solutions for T$> and T^P and thus 
022 ^ an< 3 af 3 T) (see interfacial traction continuity equations (26) and (27)), are obtained in terms of the 
inverse of the ABB'D matrix elements appearing on the left hand side of eqn (30), denoted by m, and ele- 
ments of the vectors on the right hand side of eqn (30). The knowledge of these average subcell stresses 
provides the solution for (j|^ as well. The resulting relations between the average subcell normal 
stresses, the macroscopic strains and plastic and thermal subcell strains are. 


On** = ape, 1 + M p7) e 22 + cpe» + rpAT + 4 >P 

a? = , + b $% 2 + c^e 3 3 + fpAT + dp 

c 33 7) = «f e„ + bf % 2 + cf % 3 + if >AT + Of 


where 


a p = ill-W -Spaf ] , *P =- +S®>bf>] 


cp=- 


+ Spcf ] 


(39) 

(40) 

(41) 


(42a) 


rp = - 



[ap + sp ip + 5p rp , 


•f— 5S^ } 


<Py> 


+ Spdp +Sp<Dp| ( 42b) 
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(43a) 



N y Nf 

a ^ =— ^ w (,V r +p,a)^o — X m (A , T +p,A' 7 +a) e a > 
a=l a=l 





N y N , 

*4 Y) = A X m <Y.a> , 4 T)= ^»l ( y.AVkx) 
ot=l a=l 

W T 

= X m (Y.«)^°° + Z m (Y,A' T +a)P2 a) (43b) 

a=l (x=l 

N Y N p 

bf ] =*X m (W r +p,a) . c(P> =/X m w 7 +p,/V T -Kx) (44a) 
0=1 0=1 

*T N> 

= X m (JV T +p,a)/^ a) + X m (A' 7 +p,A' T +a)P^ t) (44b) 

a=l a=l 


Substituting eqns (39)-(41) in the definitions for the average subcell normal stresses given by the 
first three of eqns (2), we obtain the macroscopic normal constitutive equations of the form given by eqn 
(6). The elements of the macroscopic elastic stiffness matrix C* , the macroscopic thermal expansion 
coefficient vector a* and the macroscopic plastic strain S P , are explicitly given below in terms of the sub- 
cell material and geometric parameters and subcell plastic strains. It is easily verified through numerical 
experiments that, despite appearance to the contrary, the stiffness matrix C* is symmetric. 


Cn Cl 2 C| 3 

C21 C22 C23 
C31 C32 C33 


1 N l 

T7X W-fiP 

nl y=l P=1 
1 N, 

7 1, 1 A* 

1 Y-l 

, N, 


hi 


M 


1 N r N » 

X X 

Y=1 P=1 

1 

yX^P 

4 Y^l 

N, 


hi 


p=l 


1 N y N f 

X XVr^i 

y=l M 
1 

yX ‘A* 

1 T-l 
N a 


1 ;v P 1 /V P 

XM* t ixvf 


^ p=l 


(45) 


* " 
<x„ 


Cn Cn C'n 

-1 

a 2 2 

► = - 

C21 C22 C23 

■< 

a 33 

V J 


Ch C32 C33 



1 

77X XWf’» 

"* -^1 p=l 

1 

yXy^ 

1 Y =1 

1 N l 

xX^rf’ 

" P=1 


r -| 




£n 


Cn Cn C*n 

-1 

£22 

y ~~ 

C21 C22 C23 

■i 

£33 

< j 


C31 C32 C33 



1 

77 X XVr 4 ^ 

ni r=i p=i 


1 


Ny 


yX^ 

1 Y=1 


1 "e 

A w— 1 


X*p^ 


P=1 


(46) 


Similarly, the relations between the average subcell shear stresses, plastic strains and macroscopic 
strains are obtained from eqns (34), (35) and (38) in a straightforward fashion in the form (taking into 
account the interfacial traction continuity conditions (32), (33) and (36), (37)), 
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( 47 ) 


Cn=-jr[hli2- P \t] 

E* y 

o ( i? = -7T-[/ei3 -p$] (48) 

F P 

^23 Y) = 77 [^^23 -P23] (49) 

Lr 

where 

£ Y =7l¥^’ M'S = 2>p£u W (50) 

z p=i p=i 

Fp = jhyS < $ ) , P® = 2^3^ (51) 

z r=i r=' 

1 Ny N> m jVy _p( By) 

G = -i-S XAp/ Y Sr , P23 = I w & i (52) 

Z Y=1 p=l T=1 P=1 


Substituting eqns (47)-(49) in the definitions for the average subcell shear stresses given by the last 
three of eqns (2), we obtain the macroscopic shear constitutive equations of the form given by eqn (6). 

* —P 

The elements of the macroscopic elastic stiffness matrix C and the macroscopic plastic strain £ are 
explicitly given in terms of the subcell material and geometric parameters, and plastic strains as follows, 


G44 0 0 

0 C55 0 
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NUMERICAL RESULTS 

We begin this section by comparing the times required by the original and reformulated versions of 
the generalized method of cells to generate the thermal heat-up response of a unidirectional silicon 
carbide/titanium aluminide (SiC/TiAl) composite. In the presence of temperature-dependent thermoelas- 
tic properties of the individual constituents, as is the case here, this is the most computationally demand- 
ing case since the solutions to the continuity equations in the original and reformulated versions, eqns 
(16) and (30) respectively, have to be generated at each thermal load increment. The temperature- 
dependent material parameters of the constituents are given in Table 1. The SiC fiber is treated as elastic 
while the titanium aluminide matrix is treated as elastic-plastic with bilinear hardening. The fiber volume 
fraction of 0.25 was employed in the calculations and the unit cell was composed of a square fiber 
embedded in a square array. The thermal heat-up response was generated for several unit cells with a pro- 
gressively greater number of subcells without changing the actual cell geometry. This was accomplished 
by subdividing each of the four subcells in the original repeating unit cell into increasingly greater 
numbers of subcells. Therefore, the thermal heat-up response obtained from repeating unit cells with dif- 
ferent subcell discretization levels was expected to remain the same for both versions of the generalized 
method of cells, which provided a check on the generated results. This is illustrated in Fig. 2. 

The comparison of the CPU times required by the two versions of the generalized method of cells to 
simulate the thermal heat up of the SiC/TiAl composite from 24°C to 815°C is provided in Table 2. These 
results were generated on an IBM RISC System/6000, Machine Type 7012 (Model 39H). The CPU times 
obtained using the original version of the generalized method of cells are given for unit cells containing 
up to 12 x 12 subcells. The CPU times for unit cells containing more than this number of subcells were 
not generated due to excessively long execution times required for their analysis. Alternatively, CPU 
times obtained using the reformulated generalized method of cells were easily determined for unit cells 
with as many as 100 x 100 subcells. In addition to the dramatic reduction in the computational time evi- 
dent in the results generated by the reformulated version, it must be pointed out that the size of the unit 
cell that can be analyzed using the original version is limited by the storage capacity of the employed 
machine. Thus while unit cells consisting of 100 x 100 subcells can readily be analyzed using the refor- 
mulated version, the maximum number of subcells employed in the original version is not nearly as large. 
This can be easily understood by recalling that while the size of the matrix that must be inverted in the 
original version is 6N pNy x 6 N^Ny, the size of the corresponding matrix in the reformulated version is 
only (N$ + Ny) x (N$ + N y ). Thus the difference in the size of the matrices that must be inverted in the 
two versions increases rapidly with increasing number of subcells. It must be mentioned that further 
reduction in the size of the strain concentration matrix in the original version could be attained by expli- 
citly taking into account the commonality of the axial strain in all the subcells, and by separating the nor- 
mal and shear contributions. However, this would not change the N$N y x N$N y order of magnitude of 
the size of the strain concentration matrix in the original version. 

We note that under certain circumstances, a further improvement in computational efficiency can be 
attained by expressing the system of equations (30) either in terms of 7^ or T^. In this case, the 
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calculation of the effective properties and thermal and plastic strains must be reformulated (the details of 
which are not included here). However, as mentioned previously, the reduction in the size of eqn (30) 
offered by this is offset by an increased number of matrix operations. For instance, for the thermal heat-up 
problem considered above, the execution times obtained using the system of equations (30) and its 
reduced counterpart are approximately the same for unit cells with up to 100 x 100 subcells. Increasing 
the number of subcells to 300 x 300 does result in an increased efficiency of the reduced system of equa- 
tions by over a factor of two (55,556 vs. 22,229 CPU seconds for the unreduced and reduced system of 
equations (30), respectively). However, this increase is obtained only for unit cells containing a suffi- 
ciently large number of subcells when the system of equations (30) is solved at every load increment, as 
in the case of temperature-dependent elastic properties of the constituent phases. When the coefficients 
appearing on the left side of eqn (30) remain constant during the loading history, as in the case of isother- 
mal mechanical loading, reducing the system of equations (30) further offers no advantage. In fact, for the 
boron/aluminum composite considered next, transverse loading of a unit cell with a square fiber contain- 
ing 300 x 300 subcells required 6,521 and 31,306 CPU seconds using the unreduced and reduced system 
of equations, respectively, almost a five-fold deterioration in the computations efficiency . 

Next, we present numerical examples that illustrate the capability of the reformulated generalized 
method of cells to efficiently model the response of unidirectional metal matrix composites with detailed 
fiber geometries. In particular, we consider circular, elliptical (with an aspect ratio of 2.5), hexagonal, 
diamond and square fibers embedded in square arrays (i.e., square repeating unit cells). To obtain true 
cross-sectional shapes of these fibers, unit cells containing 120 x 120, 68 x 68, 102 xl02, 102 x 102, and 
4x4 subcells were constructed for the circular, elliptical, hexagonal, diamond, and square fibers, respec- 
tively, Fig. 3. The properties of the fiber and matrix phases were those of an elastic boron fiber and an 
elastic-plastic aluminum matrix with bilinear hardening given in Table 3. Two hardening slopes for the 
aluminum matrix were employed in the calculations, whose magnitudes were 1/10,000 and 1/5 of the 
elastic Young’s modulus. The smaller hardening slope produces a response which is nearly elastic- 
perfectly plastic, with practically no apparent strain hardening. Figure 4 presents the stress-strain curves 
of the constituent phases employed in the calculations to generate the composite stress-strain response of 
the different repeating unit cells under transverse normal loading. Transverse normal loading was chosen 
because previous investigations have revealed that differences in the response of composites with dif- 
ferent fiber shapes are most dramatic under this manner of loading (Arnold, et al., 1996). 

Figure 5 presents the transverse response of unit cells containing 0.25 volume fraction of the five 
differently-shaped fibers embedded in the aluminum matrix with the low and the high rates of strain har- 
dening. In the case of the hexagonal and elliptical fibers, the transverse loading was applied in both of the 
principal material directions x 2 and x$ due to the cross-sectional differences along these directions. When 
the aluminum matrix exhibits no hardening, Fig. 5a, the effect of the fiber cross-section shape on the 
transverse response is only evident, and to a very limited extent, in the elastic and the initial yielding 
regions, with virtually no differences observed when the matrix is fully yielded. In the elastic region, the 
transverse response of the unit cell with the elliptical fiber is somewhat stiffer in the x 2 direction than the 
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response of the remaining unit cells, however once yielding initiates this difference vanishes. The stiffer 
response is expected since the major axis of the elliptical fiber lies along the xj direction. Alternatively, 
when the aluminum matrix strain hardens. Fig. 5b, the fiber shape has a more pronounced effect on the 
transverse response in the plastic region. Substantial differences are now observed in the and x 3 direc- 
tion responses of the unit cell with the elliptical fiber. The response of this unit cell in both directions is 
stiffer than the responses of the remaining unit cells. Among the remaining unit cells, the one with the 
square fiber exhibits the stiffest response which is slightly higher than the rest. No discernible differences 
are observed among the responses of the unit cells with diamond, hexagonal and circular fibers. 

Increasing the fiber volume fraction further accentuates the differences in the composite’s 
transverse response due to the fiber’s cross-sectional shape. Figure 6 shows the results corresponding to 
those provided in Fig. 5 when the fiber volume fraction is increased to 0.30. This fiber volume fraction is 
close to the maximum allowable for the unit cell with the elliptical fiber, which is limited by the contact 
of fibers along the major axis in two adjacent cells. This contact occurs when the fiber volume fraction is 
0.31 in the case of fibers with an aspect ratio of 2.5, which is the present situation. When the aluminum 
matrix exhibits no hardening. Fig. 6a, the effect of the fiber cross-sectional shape on the transverse 
response in the plastic region becomes discernible, with the square and elliptical fibers being the most 
effective in increasing the flow stress of the composite, and the diamond, hexagonal and circular fibers 
the least effective. When the aluminum matrix strain hardens, Fig. 6b, the nearly touching elliptical fibers 
produce the stiffest overall transverse response, with a substantially greater difference in the x 2 and x 3 
directions than that observed at the lower fiber volume fraction. Fig. 5b. The responses of the unit cells 
with the remaining fibers are substantially lower, and the differences due to the fiber cross-sectional 
shape somewhat more pronounced than in the 0.25 fiber volume fraction cases. Among these responses, 
the stiffest response is observed for the unit cell with the square fiber and the most compliant response for 
the unit cell with the circular fiber. 

Figure 7 presents the results that correspond to those shown in the preceding two figures when the 
fiber volume fraction is further increased to 0.40. This fiber volume fraction exceeds the maximum allow- 
able volume fraction for the employed elliptical fibers, and thus the results are limited to unit cells with 
square, diamond, hexagonal and circular fibers. In the case of the aluminum matrix with no hardening. 
Fig. 7a, a substantial difference between the unit cell with the square fiber and the remaining unit cells is 
now apparent in the plastic region. The square fiber provides a 50% increase in the transverse flow stress 
of the composite relative to that of the diamond and circular fibers. The transverse response of the unit 
cell with the hexagonal fiber is essentially the same for both directions of loading and lies above the 
response of the unit cells with the diamond and circular fibers. It is characterized by slight strain harden- 
ing which produces a transverse stress at the strain of 1 % that is approximately 1 8% higher than the flow 
stress of the unit cells with the diamond and circular fibers. When the aluminum matrix strain hardens, 
Fig. 7b, increasing the fiber volume fraction from 0.30 to 0.40 also increases the differences in the plastic 
regime due to the fiber cross-sectional shape. As in the preceding cases, the square fiber produces the stif- 
fest response, while the circular fiber the most compliant. The response of the unit cell with the diamond 
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fiber falls below that of the unit cell with the square fiber, while the response of the unit cell with the hex- 
agonal fiber lies just above that of the unit cell with the circular fiber for both loading directions. 

The effects of fiber shape on the transverse response of the unidirectional boron/aluminum compo- 
site illustrated herein are consistent with those reported by previous investigators that were recently sum- 
marized by Arnold, et al. (1996). In particular, the effect of fiber shape at low fiber volume fractions is 
small and increases with increasing fiber content. Square fibers provide the stiffest overall response 
(when the matrix strain hardens) and the highest flow stress (when the matrix does not strain harden) 
while circular fibers provide the most compliant response and the lowest flow stress. This is due to the 
square fiber’s ability to provide a higher magnitude of hydrostatic stress in the matrix phase relative to the 
circular fiber, thereby delaying localized yielding and providing constraint on the expansion of the plastic 
zone throughout the matrix phase. Hexagonal fibers generate an intermediate response. It must be men- 
tioned that the results reported in the literature addressing the effect of fiber shape on the response of 
MMCs are typically limited in scope, with a limited number of fiber shapes and typically a single strain 
hardening rate for the matrix phase. This is due to the computationally intensive finite-element approach 
typically employed to investigate these effects. Alternatively, the results generated with the reformulated 
version of the generalized method of cells that have been presented herein provide for the first time an 
exhaustive comparison of the fiber shape impact on the MMCs’ transverse response as a function of the 
fiber volume fraction and the matrix phase’s strain hardening capability. Similarly, the influence of the 
fiber arrangement on the response of MMCs, which has a reportedly greater impact than the fiber shape, 
can now be efficiently investigated. 

The last set of results addresses the important question of the extent to which a given fiber shape 
should be approximated (the discretization level of the repeating unit cell) in order to obtain sufficiently 
accurate results for that shape. In particular, we consider increasing levels of cell discretization for a cir- 
cular fiber to determine the level of microstructural refinement at which convergence to the "true" solu- 
tion obtained from a highly discretized unit cell is achieved. The "true" solution is based on the previ- 
ously employed unit cell containing a highly refined circular fiber shape generated using 120 x 120 sub- 
cells. The unit cells employed in this convergence study are shown in Fig. 8. In addition to the unit cells 
with the square fiber and the highly refined circular fiber, unit cells with a fiber in the shape of a cross 
containing 6x6 subcells and circular fibers approximated by 8 x 8 and 18x18 subcells were considered. 
The results for the transverse response of these unit cells containing the previously employed boron fiber 
and aluminum matrix are presented in Fig. 9 based on a fiber volume fraction of 0.40 and no strain har- 
dening of the aluminum matrix. As observed previously, the response of the unit cell containing the 
highly refined circular fiber cannot be approximated in the plastic region by that containing the square 
fiber at this fiber volume fraction. Despite the fact that the elastic responses are virtually the same, the 
substantial strain hardening produced by the square fiber leads to the macroscopic flow stress for the 
composite which is more than 50% greater than that due to the circular fiber. The cross-shaped fiber pro- 
duces a flow stress which is 20% higher than that due to the highly refined circular fiber. This difference 
reduces to 10% and 5% for the circular fibers approximated by 8 x 8 and 18x18 subcells, respectively. 
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CONCLUSIONS AND FUTURE PERSPECTIVES 

The efficient implementation of the generalized method of cells outlined herein greatly facilitates 
analysis of unidirectional composites whose microstructures require a highly detailed geometric represen- 
tation. This was accomplished by first reformulating the displacement continuity equations in the original 
formulation of the generalized method of cells in terms of subcell stress components and then directly 
incorporating the traction continuity conditions into these equations. This effectively replaces the subcell 
strains by the subcell interfacial tractions as the basic unknown microvariables, which in turn substan- 
tially reduces the size of the system of equations for the determination of these micro variables. This sub- 
stantial reduction in the size of the system of equations for the unknown microvariables makes possible 
the inelastic analysis of repeating unit cells containing thousands of subcells which could not be analyzed 
previously. Consequently, repeating unit cells with highly detailed fiber shapes can now be analyzed effi- 
ciently in the presence of matrix plasticity as demonstrated herein. Comparison of CPU times required by 
the original and reformulated versions of the generalized method of cells to simulate the thermal heat-up 
response of a SiC/TiAl composite with temperature-dependent thermoelastoplastic constituent properties 
demonstrated a dramatic reduction in the execution times as a function of the number of subcells in the 
repeating unit cell due to the reformulation. Further, while the thermal analysis of repeating unit cells 
containing more than 100 subcells proved to be impractical using the original formulation, such analysis 
was perfomed efficiently using the new formulation in a fraction of the time required of the original for- 
mulation. This opens up the possibility of using the new formulation in conjunction with large-scale 
structural analyses to investigate the response of composite structures with microstructures that require 
detailed geometrical representations. 

The examples presented herein were limited to the analysis of the transverse response of unidirec- 
tional metal matrix composites with highly refined elliptical, circular, diamond, and hexagonal fibers. 
The effect of the fiber cross section on the transverse response predicted by the reformulated generalized 
method of cells was consistent with the results reported in the literature. An additional study into the 
effect of different levels of circular fiber shape approximation on the transverse response of a undirec- 
tional composite revealed that a unit cell discretized into 18x18 subcells closely approximated the 
response obtained with a 120 x 120 subcell unit cell. The utility of the reformulated generalized method 
of cells, however, should become even more apparent when investigating the response of unidirectional 
composites with more than one type of fiber and with different fiber arrangements or arrays. The impact 
of the fiber array architecture on composite response is known to be substantially greater than that of the 
fiber shape. However, the impact of multi-phase fibers on composite response remains to be fully charac- 
terized. This can now easily be performed using the reformulated generalized method of cells. In addition, 
modeling of mechanisms that require detailed cell discretization, such as progressive damage accummula- 
tion due to crack growth or progressive interfacial debonding can also be efficiently realized. 

The outlined reformulation of the generalized method of cells in terms of subcell interfacial trac- 
tions is expected to result in an even greater computational capability enhancement when applied to the 
three-dimensional version of the generalized method of cells developed by Aboudi (1995) for modeling 
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the response of inelastic multi-phase, short-fiber composites. This version has been employed by Herako- 
vich and Baxter (1997) to investigate the effect of pore geometry on the inelastic response of porous 
media. The three-dimensional generalized method of cells has also been employed to model the mechani- 
cal response of woven graphite/copper composites by Bednarcyk, et al. (1997) using an embedded 
micromechanics- within-micromechanics model approach. In this approach, the instantaneous response of 
the individual fiber yams, which occupy specified subcells in the three-dimensional representation of the 
woven composite, is determined at each point along the loading path using the original method of cells 
(Aboudi, 1991) and then subsequently used in the global three-dimensional generalized method of cells 
analysis. However, the three-dimensional discretization of the repeating unit cell substantially increases 
the number of unknown microvariables that must be determined, thereby placing a limit on the level of 
microstructural refinement of the weave geometry. In order to be able to model three-dimensional woven 
composites using a sufficiently detailed representation of the weave geometry, the three-dimensional gen- 
eralized method of cells has recently been reformulated along similar lines presented herein for the two- 
dimensional version. The details of this reformulation and the comparison of the analytical predictions 
with experimental data will be reported in the future. 

Finally, the reformulation of the displacement continuity equations in terms of the interfacial sub- 
cell tractions as the unknown microvariables produces a system of equations characterized by a matrix 
that consists of two diagonal submatrices and two fully-populated off-diagonal submatrices in the case of 
the two-dimensional generalized method of cells. In the case of the three-dimensional generalized method 
of cells, on the other hand, the matrix consists of thirty-six submatrices of which only twelve are fully- 
populated. Therefore, further computational efficiency enhancement can be accomplished for both ver- 
sions of the generalized method of cells through Orozco’s sparse implementation approach. 
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APPENDIX 


Relations between subcell normal tractions and macroscopic normal strains 

To develop the relations between subcell normal tractions and macroscopic normal strains, we first 
express the subcell normal strains in terms of the subcell normal stresses, and plastic and thermal strains 
using eqn (25). 
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Eauations (A5) and (A6) are subsequently used in the interfacial displacement continuity conditions (12b) 
~ — (Py) , — (Py) 

and (12d) in order to express them in terms of the normal subcell stresses o 2 2 and o 33 , 
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Use of the interfacial traction continuity conditions (26) and (27) in the above expresions yields eqns (28) 
and (29). 

Relations between subcell shear tractions and macroscopic shear strains 

To develop the relations between subcell shear tractions and macroscopic shear strains, we first 
express the subcell shear strains in terms of the subcell shear stresses, and plastic and thermal strains 
using eqn (25). 
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Equations (A9) through (All) are subsequently used in the interfacial displacement continuity conditions 
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Use of the interfacial traction continuity conditions (32), (33), (36) and (37) in the above expresions 
yields eqns (34), (35) and (38). 
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Table 1 . Material properties of SiC fiber and titanium aluminide matrix. 


Material properties 

24 °C 

200 °C 

425 °C 

600 °C 

650 °C 

815 °C 

SiC fiber 

a (x 1CT 6 cm/cm/ °C) 

3.53 

3.62 

3.87 

4.19 

4.28 

4.5 

Young’s modulus (GPa) 

400.0 

400.0 

400.0 

400.0 

400.0 

400.0 

Poisson’s ratio 

0.25 

0.25 

0.25 

0.25 

0.25 

0.25 

Ti-24A1- 1 1 Nb matrix 

a (x 10 -6 cm/cm/ °C) 

9.0 

9.36 

10.26 

10.53 

10.62 

11.07 

Young’s modulus (GPa) 

110.3 

100.0 

75.8 

86.2 

68.2 

11.2 

Poisson’s ratio 

0.26 

0.26 

0.26 

0.26 

0.26 

0.26 

Yield stress (MPa) 

371.5 

406.7 

370.2 

290.9 

269.5 

165.5 

Hardening slope (GPa) 

22.98 

3.04 

2.22 

1.29 

0.67 

0.00 


Table 2. Comparison of CPU times required by the original and reformulated versions of the generalized 
method of cells to generate a thermal heat-up response of a SiC/TiAl unidirectional composite. 


GMC version CPU times (seconds) as a function of the number of subcells 



2x2 

4x4 

6x6 

8x8 

10x10 

12x12 

20x20 

100x100 

original 

0.87 

19 

182 

508 

8,679 

43,781 

- 

- 

reformulated 

0.18 

0.25 

0.5 

0.9 

1.5 

2.3 

8.3 

796 


Table 3. Material properties of the boron fiber and aluminum matrix. 


Material 

Young’s modulus (GPa) 

Poisson’s ratio 

Hardening slope (GPa) 

Boron fiber 

414.0 

0.20 

N/A 

Aluminum matrix 

55.2 

0.30 

E/ 10,000 and E/5 


Note: The hardening slope is the secondary slope in a bilinear stress-strain representation of the matrix 
phase’s elastic-plastic response. 
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x 2 



Doubly periodic array of fibers 



Discretization of the repeating unit cell 


Figure 1 . A unidirectional composite modeled by the generalized method of cells as a doubly-periodic 
array of fibers (top), and the details of the repeating unit cell discretization (bottom). 
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Figure 2. Thermal response of a unidirectional SiC/TiAl composite predicted by the original and 
reformulated versions of the generalized method of cells: (a) longitudinal; (b) transverse. 
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Figure 5. Transverse response of a unidirectional B/Al composite with differently-shaped fibers 
containing 0.25 fiber volume fraction: (a) H = E/10,000; (b) H = E/5. 
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Figure 7. Transverse response of a unidirectional B/Al composite with differently-shaped fibers 
containing 0.40 fiber volume fraction: (a) H = E/1 0,000; (b) H = E/5. 



Figure 8. Repeating unit cells containing fibers with increasingly refined circular fiber shape. 


32 






Transverse Normal Stress (Ml 



Strain (%) 


Figure 9. Transverse response of repeating unit cells containing fibers with increasingly refined circular 

fiber shape. 
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